Imports System.Windows.Forms
Imports Crownwood.Magic.Docking
Imports Crownwood.Magic.Controls
Imports Crownwood.Magic.Common
Imports Crownwood.Magic.Menus

Public Class SegmentationForm
    'Dim OriginalImage As Bitmap

    Sub New()

        ' This call is required by the Windows Form Designer.

        InitializeComponent()

    End Sub

    Sub New(ByVal Bit As Bitmap, ByVal max As Integer, ByVal min As Integer)

        '   OriginalImage = B

        ' This call is required by the Windows Form Designer.


        
        InitializeComponent()
        PictureBox3.Image = Bit
        OriginalImage = Bit
        ' Add any initialization after the InitializeComponent() call.
        ValueSlider.Min = min
        ValueSlider.Max = max
        MinTextBox.Text = CStr(65)
        MaxTextBox.Text = CStr(128)

        ' Add any initialization after the InitializeComponent() call.

        Me.ShowIcon = False

    End Sub

    Dim SegmentationClassObject As New ExpectationMaximization

    Private Sub okbutton_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles OKButton.Click

        Me.DialogResult = System.Windows.Forms.DialogResult.OK


        Dim CSF_GM, GM_WM As Int16
        CSF_GM = Convert.ToInt16((ValueSlider.Min / 255) * (FormMain.AnatomyVolume.Matrix.Maximum - FormMain.AnatomyVolume.Matrix.Minimum) + FormMain.AnatomyVolume.Matrix.Minimum)
        GM_WM = Convert.ToInt16((ValueSlider.Max / 255) * (FormMain.AnatomyVolume.Matrix.Maximum - FormMain.AnatomyVolume.Matrix.Minimum) + FormMain.AnatomyVolume.Matrix.Minimum)

        SegmentationClassObject = New ExpectationMaximization(FormMain.AnatomyVolume.Matrix, CSF_GM, GM_WM)
        Dim count As Integer = 1
        While SegmentationClassObject.ComputePosteriorProbabilitiesAndUpdate(True) = False
            count += 1
            SegmentationClassObject.ComputePosteriorProbabilitiesAndUpdate(True)
        End While

        SegmentationClassObject.DetectEdges()
        MsgBox("Done. Iterations : " & count & ".")
        CheckedListBox1.Enabled = True
        GroupBox1.Visible = True
        PictureBox3.Image = SegmentationClassObject.ClassifierValues.WM.SliceToBitmap(SegmentationClassObject.ClassifierValues.WM.TransverseSlice(CInt(SegmentationClassObject.ClassifierValues.GM.z / 2)))

    End Sub

    Sub CurrentSliceInformation(ByVal sender As Object, ByVal e As EventArgs)
    End Sub

    Private Sub mintextbox_TextChanged(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MinTextBox.TextChanged
        ValueSlider.Min = Byte.Parse(MinTextBox.Text)
        refreshpreview()
    End Sub

    Private Sub maxtextbox_TextChanged(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MaxTextBox.TextChanged
        ValueSlider.Max = Byte.Parse(MaxTextBox.Text)
        refreshpreview()
    End Sub

    Private Sub ValueSlider_ValuesChanged(ByVal sender As Object, ByVal e As System.EventArgs) Handles ValueSlider.ValuesChanged
        MinTextBox.Text = ValueSlider.Min.ToString
        MaxTextBox.Text = ValueSlider.Max.ToString
    End Sub

    Sub refreshpreview()

        Dim temp As New Bitmap(FormMain.AnatomyDisplay.CurrentImage(False))
        For i As Integer = 0 To FormMain.AnatomyDisplay.CurrentImage(False).Width - 1
            For j As Integer = 0 To FormMain.AnatomyDisplay.CurrentImage(False).Height - 1
                If temp.GetPixel(i, j).R < ValueSlider.Min Then
                    temp.SetPixel(i, j, Color.FromArgb(0, 0, 0)) 'noise and csf
                ElseIf temp.GetPixel(i, j).R > ValueSlider.Max Then
                    temp.SetPixel(i, j, Color.FromArgb(255, 255, 255)) 'wm
                Else
                    temp.SetPixel(i, j, Color.FromArgb(128, 128, 128)) 'gm
                End If
            Next
        Next
        PictureBox3.Image = temp
        modifiedImage = temp
    End Sub

    Dim modifiedImage, OriginalImage As Bitmap
    Private Sub picturebox3_MouseDown(ByVal sender As Object, ByVal e As System.Windows.Forms.MouseEventArgs) Handles PictureBox3.MouseDown
        modifiedImage = CType(PictureBox3.Image, Bitmap)
        PictureBox3.Image = OriginalImage
        PictureBox3.Refresh()
    End Sub

    Private Sub picturebox3_MouseUp(ByVal sender As Object, ByVal e As System.Windows.Forms.MouseEventArgs) Handles PictureBox3.MouseUp
        PictureBox3.Image = modifiedImage
        PictureBox3.Refresh()
    End Sub

    Private Sub CloseButton_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles CloseButton.Click
        Me.Close()
    End Sub

   
    'Private Sub Button7_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button7.Click
    '    'laplace smoothen
    '    Dim a As InputBoxResult = CustomInputBox.Show("Iterations, Corner Cutoff Fraction:" & _
    '    vbNewLine & vbNewLine & "Note: Iterations defines the number of times" & _
    '                vbNewLine & "the laplace smoothing is applied. The Cutoff" & _
    '                vbNewLine & "defines the minimum zeros a voxel should have" & _
    '                vbNewLine & "around it to avoid getting smoothened.", "Parameter", "100,0.5")
    '    If a.ReturnCode = Windows.Forms.DialogResult.OK Or a.ReturnCode = Windows.Forms.DialogResult.Yes Then
    '        Dim temp, count, zero As Long
    '        For c As Integer = 0 To CInt(CDbl(a.Text.Split(CChar(","))(0)) - 1)
    '            Me.Text = "Segmentation :: Cortical Smoothing " & " - " & CInt(c * 100 / CInt(CDbl(a.Text.Split(CChar(","))(0)) - 1)) & "% Complete."
    '            For i As Integer = 1 To GM.x - 2
    '                For j As Integer = 1 To GM.y - 2
    '                    For k As Integer = 1 To GM.z - 2
    '                        If GM(i, j, k) <> 0 Then 'skip 0 and boundary
    '                            temp = 0
    '                            count = 0
    '                            zero = 0
    '                            For ki As Integer = -1 To 1 Step 1
    '                                For kj As Integer = -1 To 1 Step 1
    '                                    For kk As Integer = -1 To 1 Step 1
    '                                        If GM(i + ki, j + kj, k + kk) <> 0 Then
    '                                            count = CShort(count + 1)
    '                                            temp += GM(i + ki, j + kj, k + kk)
    '                                        Else
    '                                            zero = CShort(zero + 1)
    '                                        End If
    '                                    Next
    '                                Next
    '                            Next
    '                            temp = CShort(temp / count)
    '                            If ((zero / 26) <= CDbl(a.Text.Split(CChar(","))(1))) Then GM(i, j, k) = CShort(temp)
    '                        End If
    '                    Next
    '                Next
    '            Next
    '            PictureBox3.Image = GM.TransverseSliceB(CInt(GM.z / 2))
    '            PictureBox3.Refresh()
    '        Next
    '        FormMain.AnatomyDisplay.OpenFile(GM, "Smoothened GM")
    '        temp = 0
    '        count = 0
    '        Dim Nx, Ny, Nz, N, Length As New MatrixDataStructures.Matrix3DSingle(GM.x, GM.y, GM.z)
    '        For i As Integer = 2 To GM.x - 2
    '            Me.Text = "Segmentation :: Computing Gradient " & " - " & (CInt(i * 100 / GM.x) - 1) & "% Complete."
    '            For j As Integer = 2 To GM.y - 2
    '                For k As Integer = 2 To GM.z - 2
    '                    If GM(i, j, k) > 0 Then
    '                        temp += 1
    '                        Nx(i, j, k) = -GM(i, j, k) + GM(i - 1, j, k)
    '                        Ny(i, j, k) = -GM(i, j, k) + GM(i, j - 1, k)
    '                        Nz(i, j, k) = -GM(i, j, k) + GM(i, j, k - 1)
    '                        N(i, j, k) = (((((Nx(i, j, k) ^ 2) + (Ny(i, j, k) ^ 2) + (Nz(i, j, k) ^ 2)) ^ 0.5)))
    '                        Nx(i, j, k) = (Nx(i, j, k) / N(i, j, k))
    '                        Ny(i, j, k) = (Ny(i, j, k) / N(i, j, k))
    '                        Nz(i, j, k) = (Nz(i, j, k) / N(i, j, k))
    '                        If GM(i, j, k) >= 248 Then count += 1
    '                    End If
    '                Next
    '            Next
    '        Next
    '        MsgBox("Total Points: " & temp & ". Total Seeds: " & count)

    '        Dim dt As Double = 0.1
    '        Dim nextx As Integer = 0
    '        Dim nexty As Integer = 0
    '        Dim nextz As Integer = 0

    '        FormMain.DiffusionDisplay.OpenFile(N, "1e")
    '        FormMain.FAVolume.MatrixD = N
    '        N.Dispose()

    '        For i As Integer = 2 To GM.x - 2
    '            Me.Text = "Segmentation :: Computing Direction Lengths " & " - " & (CInt(i * 100 / GM.x) - 1) & "% Complete."
    '            For j As Integer = 2 To GM.y - 2
    '                For k As Integer = 2 To GM.z - 2
    '                    Length(i, j, k) = 0
    '                    If GM(i, j, k) >= 240 Then
    '                        'Euler's Method
    '                        nextx = CInt(i + Nx(i, j, k) * dt)
    '                        nexty = CInt(j + Ny(i, j, k) * dt)
    '                        nextz = CInt(k + Nz(i, j, k) * dt)
    '                        While GM(nextx, nexty, nextz) > 0 And GM(nextx, nexty, nextz) <> 250
    '                            Length(i, j, k) += 1
    '                            nextx = CInt(nextx + Nx(i, j, k) * dt)
    '                            nexty = CInt(nexty + Ny(i, j, k) * dt)
    '                            nextz = CInt(nextz + Nz(i, j, k) * dt)
    '                            If nextx < 0 Or nextx > GM.x - 1 Then Exit While
    '                            If nexty < 0 Or nexty > GM.y - 1 Then Exit While
    '                            If nextz < 0 Or nextz > GM.z - 1 Then Exit While
    '                        End While
    '                        nextx = CInt(i - Nx(i, j, k) * dt)
    '                        nexty = CInt(j - Ny(i, j, k) * dt)
    '                        nextz = CInt(k - Nz(i, j, k) * dt)
    '                        While GM(nextx, nexty, nextz) > 0 And GM(nextx, nexty, nextz) <> 250
    '                            Length(i, j, k) += 1
    '                            nextx = CInt(nextx - Nx(i, j, k) * dt)
    '                            nexty = CInt(nexty - Ny(i, j, k) * dt)
    '                            nextz = CInt(nextz - Nz(i, j, k) * dt)
    '                            If nextx < 0 Or nextx > GM.x - 1 Then Exit While
    '                            If nexty < 0 Or nexty > GM.y - 1 Then Exit While
    '                            If nextz < 0 Or nextz > GM.z - 1 Then Exit While
    '                        End While
    '                    End If
    '                Next
    '            Next
    '        Next
    '        MsgBox("!")

    '        FormMain.AnatomyDisplay.OpenFile(Length, "Length")
    '    End If
    'End Sub

    Private Sub CheckedListBox1_SelectedIndexChanged(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles CheckedListBox1.SelectedIndexChanged
        Select Case CheckedListBox1.SelectedIndex
            Case 0
                FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(SegmentationClassObject.ReturnWMData, FormMain.AnatomyVolume), "WM", True, False)
            Case 1
                FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(SegmentationClassObject.ReturnGMData, FormMain.AnatomyVolume), "GM", True, False)
            Case 2
                FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(SegmentationClassObject.ReturnCSFData, FormMain.AnatomyVolume), "CSF", True, False)
            Case 3
                FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(SegmentationClassObject.ReturnWMMap, FormMain.AnatomyVolume), "WM", True, False)
            Case 4
                FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(SegmentationClassObject.ReturnGMMap, FormMain.AnatomyVolume), "GM", True, False)
            Case 5
                FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(SegmentationClassObject.ReturnCSFMap, FormMain.AnatomyVolume), "CSF", True, False)
            Case 6
                SegmentationClassObject.PromptSaveWM()
            Case 7
                SegmentationClassObject.PromptSaveGM()
            Case 8
                SegmentationClassObject.PromptSaveCSF()
            Case 9
                FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(Me.SegmentationClassObject.ClassifierValues.WMGMEdge, FormMain.AnatomyVolume), "GM-WM Edge", True, False)
            Case 10
                FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(SegmentationClassObject.ClassifierValues.GMCSFEdge, FormMain.AnatomyVolume), "GM-CSF Edge", True, False)
        End Select
        PictureBox3.Image = FormMain.AnatomyDisplay.CurrentImage(False)
        PictureBox3.Refresh()

    End Sub

   
    Private Sub Button7_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button7.Click
        SegmentationClassObject.ApplyLaplaceSmoothing(100)
        PictureBox3.Image = FormMain.AnatomyDisplay.CurrentImage(False)
        PictureBox3.Refresh()
    End Sub

    Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click
        SegmentationClassObject.ComputeGradients()
    End Sub
End Class

Public Class ExpectationMaximization
    Public Structure ClassificationInformation
        Public GMMean, WMMean, CSFMean, GMVar, WMVar, CSFVar, GMPrior, WMPrior, CSFPrior As Double
        Public WM, GM, CSF, Data, Matrix, WMGMEdge, GMCSFEdge As MatrixDataStructures.Matrix3DInt16
        Public countGM, countWM, countCSF As Long

    End Structure
    Dim History As String
    Public ClassifierValues As ClassificationInformation
    Public LaplaceField As MatrixDataStructures.Matrix3DInt16
    ' Steps for EM Based Classification
    '-----------------------------------
    ' Compute initial Mean, Var & Prior Probability of GM WM CSF from User's threshold
    ' Model each class as gaussian centered at the class mean
    ' Compute normalized posterior probability (Gaussian * prior)
    ' Compute new Var & Prior Probabilities and Classify Again
    ' Keep going till new Var & prior = old 

    Public Sub New()
        ' default constructor
    End Sub

    Public Sub New(ByVal UnSegmentedData As MatrixDataStructures.Matrix3DInt16, ByVal GM_CSF_Threshold As Integer, ByVal WM_GM_Threshold As Integer)

        ClassifierValues.WM = New MatrixDataStructures.Matrix3DInt16(UnSegmentedData.x, UnSegmentedData.y, UnSegmentedData.z)
        ClassifierValues.GM = New MatrixDataStructures.Matrix3DInt16(UnSegmentedData.x, UnSegmentedData.y, UnSegmentedData.z)
        ClassifierValues.CSF = New MatrixDataStructures.Matrix3DInt16(UnSegmentedData.x, UnSegmentedData.y, UnSegmentedData.z)
        ClassifierValues.Data = New MatrixDataStructures.Matrix3DInt16(UnSegmentedData.x, UnSegmentedData.y, UnSegmentedData.z)
        'MsgBox(GM_CSF_Threshold & "," & WM_GM_Threshold)
        Dim sumGM As ULong = 0
        Dim sumWM As ULong = 0
        Dim sumCSF As ULong = 0
        Dim sumGM2 As ULong = 0
        Dim sumWM2 As ULong = 0
        Dim sumCSF2 As ULong = 0
        For i As Integer = 0 To UnSegmentedData.x - 1
            For j As Integer = 0 To UnSegmentedData.y - 1
                For k As Integer = 0 To UnSegmentedData.z - 1
                    If UnSegmentedData(i, j, k) > GM_CSF_Threshold / 10 Then ClassifierValues.Data(i, j, k) = UnSegmentedData(i, j, k)
                    If ClassifierValues.Data(i, j, k) > GM_CSF_Threshold / 5 And ClassifierValues.Data(i, j, k) < WM_GM_Threshold * 2 Then
                        If ClassifierValues.Data(i, j, k) <= GM_CSF_Threshold Then
                            ClassifierValues.CSF(i, j, k) = 1
                            sumCSF = CULng(sumCSF + ClassifierValues.Data(i, j, k))
                            sumCSF2 = CULng(sumCSF2 + (ClassifierValues.Data(i, j, k) * ClassifierValues.Data(i, j, k)))
                            ClassifierValues.countCSF += 1
                        ElseIf ClassifierValues.Data(i, j, k) >= WM_GM_Threshold Then
                            ClassifierValues.WM(i, j, k) = 1
                            sumWM = CULng(sumWM + ClassifierValues.Data(i, j, k))
                            sumWM2 = CULng(sumWM2 + ClassifierValues.Data(i, j, k) * ClassifierValues.Data(i, j, k))
                            ClassifierValues.countWM += 1
                        Else
                            ClassifierValues.GM(i, j, k) = 1
                            sumGM = (sumGM + CULng(ClassifierValues.Data(i, j, k)))
                            sumGM2 = sumGM2 + CULng(ClassifierValues.Data(i, j, k) * ClassifierValues.Data(i, j, k))
                            ClassifierValues.countGM += 1
                        End If
                    End If
                Next
            Next
        Next

        ClassifierValues.CSFMean = (sumCSF / ClassifierValues.countCSF)
        ClassifierValues.WMMean = (sumWM / ClassifierValues.countWM)
        ClassifierValues.GMMean = (sumGM / ClassifierValues.countGM)


        ClassifierValues.GMVar = ((sumGM2 - (sumGM * ClassifierValues.GMMean)) / (ClassifierValues.countGM - 1))
        ClassifierValues.WMVar = ((sumWM2 - (sumWM * ClassifierValues.WMMean)) / (ClassifierValues.countWM - 1))
        ClassifierValues.CSFVar = ((sumCSF2 - (sumCSF * ClassifierValues.CSFMean)) / (ClassifierValues.countCSF - 1))

        ClassifierValues.GMPrior = (ClassifierValues.countGM / (ClassifierValues.countGM + ClassifierValues.countWM + ClassifierValues.countCSF))
        ClassifierValues.WMPrior = (ClassifierValues.countWM / (ClassifierValues.countGM + ClassifierValues.countWM + ClassifierValues.countCSF))
        ClassifierValues.CSFPrior = (ClassifierValues.countCSF / (ClassifierValues.countGM + ClassifierValues.countWM + ClassifierValues.countCSF))

        History = ClassifierValues.GMMean & ", Prior" & ClassifierValues.GMPrior & ",variance (GM)" & ClassifierValues.GMVar & vbNewLine & _
        ClassifierValues.WMMean & "," & ClassifierValues.WMPrior & ",(WM)" & ClassifierValues.WMVar & vbNewLine & ClassifierValues.CSFMean & "," & ClassifierValues.CSFPrior & ",(CSF)" & ClassifierValues.CSFVar

    End Sub

    Public Function ComputePosteriorProbabilitiesAndUpdate(Optional ByVal UpdateMean As Boolean = False) As Boolean
        Dim ans As Boolean = True
        Dim a, b, c As Double
        Dim amean, bmean, cmean, avar, bvar, cvar, ameancount, bmeancount, cmeancount As Double
        For i As Integer = 0 To ClassifierValues.Data.x - 1
            For j As Integer = 0 To ClassifierValues.Data.y - 1
                For k As Integer = 0 To ClassifierValues.Data.z - 1

                    ClassifierValues.WM(i, j, k) = 0
                    ClassifierValues.GM(i, j, k) = 0
                    ClassifierValues.CSF(i, j, k) = 0

                    If ClassifierValues.Data(i, j, k) > 0 Then
                        a = Gaussian(ClassifierValues.Data(i, j, k), ClassifierValues.WMMean, ClassifierValues.WMVar) * ClassifierValues.WMPrior
                        b = Gaussian(ClassifierValues.Data(i, j, k), ClassifierValues.GMMean, ClassifierValues.GMVar) * ClassifierValues.GMPrior
                        c = Gaussian(ClassifierValues.Data(i, j, k), ClassifierValues.CSFMean, ClassifierValues.CSFVar) * ClassifierValues.CSFPrior


                        If a >= Math.Max(b, c) Then
                            amean = (amean + ClassifierValues.Data(i, j, k))
                            avar = (avar + (ClassifierValues.Data(i, j, k) * ClassifierValues.Data(i, j, k)))
                            ameancount += 1
                            ClassifierValues.WM(i, j, k) = 1
                        ElseIf b >= Math.Max(a, c) Then
                            bmean = (bmean + ClassifierValues.Data(i, j, k))
                            bvar = (bvar + ClassifierValues.Data(i, j, k) * ClassifierValues.Data(i, j, k))
                            bmeancount += 1
                            ClassifierValues.GM(i, j, k) = 1
                        ElseIf c >= Math.Max(a, b) Then
                            cmean = (cmean + ClassifierValues.Data(i, j, k))
                            cvar = (cvar + ClassifierValues.Data(i, j, k) * ClassifierValues.Data(i, j, k))
                            cmeancount += 1
                            ClassifierValues.CSF(i, j, k) = 1
                        End If
                    End If
                Next
            Next
        Next


        ClassifierValues.CSFPrior = (cmeancount / (ameancount + bmeancount + cmeancount))
        ClassifierValues.WMPrior = (ameancount / (ameancount + bmeancount + cmeancount))
        ClassifierValues.GMPrior = (bmeancount / (ameancount + bmeancount + cmeancount))
        
        If UpdateMean Then
            'compute mean 
            ClassifierValues.WMMean = (amean / ameancount)
            ClassifierValues.GMMean = (bmean / bmeancount)
            ClassifierValues.CSFMean = (cmean / cmeancount)
        End If

        ClassifierValues.WMVar = Math.Abs((avar - (ClassifierValues.WMMean * amean)) / (ameancount - 1))
        ClassifierValues.GMVar = Math.Abs((bvar - (ClassifierValues.GMMean * bmean)) / (bmeancount - 1))
        ClassifierValues.CSFVar = Math.Abs((cvar - (ClassifierValues.CSFMean * cmean)) / (cmeancount - 1))

        Dim historynew As String = ClassifierValues.GMMean & ", Prior" & ClassifierValues.GMPrior & ",variance (GM)" & ClassifierValues.GMVar & vbNewLine & _
        ClassifierValues.WMMean & "," & ClassifierValues.WMPrior & ",(WM)" & ClassifierValues.WMVar & vbNewLine & ClassifierValues.CSFMean & "," & ClassifierValues.CSFPrior & ",(CSF)" & ClassifierValues.CSFVar

        If historynew = History Then
            Dim o As New MorphologyFilters
            'Me.ClassifierValues.GM = o.AreaOpen(Me.ClassifierValues.GM, 10)
            'Me.ClassifierValues.GM = o.AreaOpen(Me.ClassifierValues.WM, 10)
            'Me.ClassifierValues.GM = o.AreaOpen(Me.ClassifierValues.CSF, 10)

            Return True
        Else
            History = historynew
            Return False
        End If
        
    End Function

    Private Function Gaussian(ByVal Value As Double, ByVal Mean As Double, ByVal Variance As Double) As Double
        Return ((1 / Math.Sqrt(2 * Math.PI * Variance)) * Math.Exp((Mean - Value) ^ 2) / (-2 * Variance))
    End Function

    Public Sub DetectEdges()
        ClassifierValues.WMGMEdge = New MatrixDataStructures.Matrix3DInt16(Me.ClassifierValues.Data.x, Me.ClassifierValues.Data.y, Me.ClassifierValues.Data.z)
        ClassifierValues.GMCSFEdge = New MatrixDataStructures.Matrix3DInt16(Me.ClassifierValues.Data.x, Me.ClassifierValues.Data.y, Me.ClassifierValues.Data.z)

        Dim GMWMEroded As New MatrixDataStructures.Matrix3DInt16(ClassifierValues.Data.x, ClassifierValues.Data.y, ClassifierValues.Data.z)

        Dim a As New MorphologyFilters
        Me.ClassifierValues.GM = a.AreaOpen(Me.ClassifierValues.GM, 10)

        Dim DilatedWM As MatrixDataStructures.Matrix3DInt16 = a.Dilate(Me.ClassifierValues.WM)
        For i As Integer = 0 To ClassifierValues.WMGMEdge.x - 1
            For j As Integer = 0 To ClassifierValues.WMGMEdge.y - 1
                For k As Integer = 0 To ClassifierValues.WMGMEdge.z - 1
                    If ClassifierValues.GM(i, j, k) > 0 Or ClassifierValues.WM(i, j, k) > 0 Then GMWMEroded(i, j, k) = 1
                    If ClassifierValues.GM(i, j, k) > 0 And DilatedWM(i, j, k) > 0 Then ClassifierValues.WMGMEdge(i, j, k) = 1
                Next
            Next
        Next
        GMWMEroded = a.Erode(GMWMEroded)
        DilatedWM.Dispose()
        For i As Integer = 0 To ClassifierValues.WMGMEdge.x - 2
            For j As Integer = 0 To ClassifierValues.WMGMEdge.y - 2
                For k As Integer = 0 To ClassifierValues.WMGMEdge.z - 2
                    If ClassifierValues.GM(i, j, k) > 0 And GMWMEroded(i, j, k) = 0 Then ClassifierValues.GMCSFEdge(i, j, k) = 1
                Next
            Next
        Next

    End Sub

    Sub ApplyLaplaceSmoothing(Optional ByVal Iterations As Integer = 1)

        Dim Sum As Long
        Dim Count As Integer
        Dim zero As Integer

        Static Data As New MatrixDataStructures.Matrix3DInt16(ClassifierValues.GM.x, ClassifierValues.GM.y, ClassifierValues.GM.z)
        For i As Integer = 1 To Data.x - 2
            For j As Integer = 1 To Data.y - 2
                For k As Integer = 1 To Data.z - 2
                    If ClassifierValues.GM(i, j, k) > 0 Then Data(i, j, k) = 10
                    If ClassifierValues.WMGMEdge(i, j, k) > 0 Then Data(i, j, k) = 240
                Next
            Next
        Next

        For c As Integer = 0 To Iterations - 1
            For i As Integer = 1 To Data.x - 2
                For j As Integer = 1 To Data.y - 2
                    For k As Integer = 1 To Data.z - 2
                        Sum = 0
                        Count = 0
                        zero = 0
                        If Data(i, j, k) > 0 And ClassifierValues.GMCSFEdge(i, j, k) = 0 And ClassifierValues.WMGMEdge(i, j, k) = 0 Then
                            If Data(i - 1, j, k) > 0 Then
                                Sum += Data(i - 1, j, k)
                                Count += 1
                            End If
                            If Data(i + 1, j, k) > 0 Then
                                Sum += Data(i + 1, j, k)
                                Count += 1
                            End If
                            If Data(i, j - 1, k) > 0 Then
                                Sum += Data(i, j - 1, k)
                                Count += 1
                            End If
                            If Data(i, j + 1, k) > 0 Then
                                Sum += Data(i, j + 1, k)
                                Count += 1
                            End If
                            If Data(i, j, k - 1) > 0 Then
                                Sum += Data(i, j, k - 1)
                                Count += 1
                            End If
                            If Data(i, j, k + 1) > 0 Then
                                Sum += Data(i, j, k + 1)
                                Count += 1
                            End If
                            If Count > 0 Then Data(i, j, k) = CShort(Sum / Count)
                        End If
                    Next
                Next
            Next
        Next
        FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(Data, FormMain.AnatomyVolume), "Laplace Vector Field", True, True)
        LaplaceField = New MatrixDataStructures.Matrix3DInt16(Data.data)
    End Sub

    Sub ComputeGradients()
        Dim GradientX, GradientY, GradientZ, GradientMax, GradientAbs As New MatrixDataStructures.Matrix3DInt16(LaplaceField.x, LaplaceField.y, LaplaceField.z)
        For i As Integer = 1 To LaplaceField.x - 2
            For j As Integer = 1 To LaplaceField.y - 2
                For k As Integer = 1 To LaplaceField.z - 2
                    If LaplaceField.data(i, j, k) > 0 Then
                        GradientX.data(i, j, k) = CShort(0.5 * (LaplaceField(i + 1, j, k) - LaplaceField(i - 1, j, k))) 'dx
                        GradientY.data(i, j, k) = CShort(0.5 * (LaplaceField(i, j + 1, k) - LaplaceField(i, j - 1, k))) 'dy
                        GradientZ.data(i, j, k) = CShort(0.5 * (LaplaceField(i, j, k + 1) - LaplaceField(i, j, k - 1))) 'dz
                        GradientAbs.data(i, j, k) = Math.Abs(GradientX.data(i, j, k)) + Math.Abs(GradientY.data(i, j, k) + Math.Abs(GradientZ.data(i, j, k)))
                        If GradientAbs.data(i, j, k) <> 0 Then 'normalize
                            GradientX.data(i, j, k) = CShort(GradientX.data(i, j, k) / ((GradientX.data(i, j, k) ^ 2 + GradientY.data(i, j, k) ^ 2 + GradientZ.data(i, j, k) ^ 2) ^ 0.5))
                            GradientY.data(i, j, k) = CShort(GradientY.data(i, j, k) / ((GradientX.data(i, j, k) ^ 2 + GradientY.data(i, j, k) ^ 2 + GradientZ.data(i, j, k) ^ 2) ^ 0.5))
                            GradientZ.data(i, j, k) = CShort(GradientZ.data(i, j, k) / ((GradientX.data(i, j, k) ^ 2 + GradientY.data(i, j, k) ^ 2 + GradientZ.data(i, j, k) ^ 2) ^ 0.5))
                            GradientAbs.data(i, j, k) = Math.Abs(GradientX.data(i, j, k)) + Math.Abs(GradientY.data(i, j, k) + Math.Abs(GradientZ.data(i, j, k)))
                            GradientMax.data(i, j, k) = Math.Max(GradientX.data(i, j, k), Math.Max(GradientY.data(i, j, k), GradientZ.data(i, j, k)))
                        End If
                    End If
                Next
            Next
        Next
        FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(GradientAbs, FormMain.AnatomyVolume), "Gradient Thickness", True, True)

        Dim OriginalVector, FinalVector As New MatrixDataStructures.Coordinate(0, 0, 0)
        Dim Thickness As New MatrixDataStructures.Matrix3DSingle(LaplaceField.x, LaplaceField.y, LaplaceField.z)
        For i As Integer = 1 To LaplaceField.x - 2
            For j As Integer = 1 To LaplaceField.y - 2
                For k As Integer = 1 To LaplaceField.z - 2
                    If Me.ClassifierValues.WMGMEdge(i, j, k) > 0 Then
                        OriginalVector.x = CShort(i)
                        OriginalVector.y = CShort(j)
                        OriginalVector.z = CShort(k)
                        Dim count As Integer = 0
                        Do While GradientAbs.data(OriginalVector.x, OriginalVector.y, OriginalVector.z) > 0 And Me.ClassifierValues.GMCSFEdge(i, j, k) = 0 And count < 20
                            Thickness(i, j, k) += 1
                            count += 1
                            FinalVector.x = OriginalVector.x + GradientX(OriginalVector.x, OriginalVector.y, OriginalVector.z)
                            FinalVector.y = OriginalVector.y + GradientY(OriginalVector.x, OriginalVector.y, OriginalVector.z)
                            FinalVector.z = OriginalVector.z + GradientZ(OriginalVector.x, OriginalVector.y, OriginalVector.z)
                            OriginalVector.x = FinalVector.x
                            OriginalVector.y = FinalVector.y
                            OriginalVector.z = FinalVector.z
                        Loop

                    End If
                Next
            Next
        Next
        FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(Thickness, FormMain.AnatomyVolume), "Step (WM-GM Edge) Thickness", True, True)
        Thickness = New MatrixDataStructures.Matrix3DSingle(LaplaceField.x, LaplaceField.y, LaplaceField.z)
        For i As Integer = 1 To LaplaceField.x - 2
            For j As Integer = 1 To LaplaceField.y - 2
                For k As Integer = 1 To LaplaceField.z - 2
                    If Me.ClassifierValues.GMCSFEdge(i, j, k) > 0 Then
                        OriginalVector.x = CShort(i)
                        OriginalVector.y = CShort(j)
                        OriginalVector.z = CShort(k)
                        Dim count As Integer
                        Do While GradientAbs.data(OriginalVector.x, OriginalVector.y, OriginalVector.z) > 0 And Me.ClassifierValues.WMGMEdge(i, j, k) = 0 And count < 20
                            count += 1
                            Thickness(i, j, k) += 1
                            FinalVector.x = OriginalVector.x + GradientX(OriginalVector.x, OriginalVector.y, OriginalVector.z)
                            FinalVector.y = OriginalVector.y + GradientY(OriginalVector.x, OriginalVector.y, OriginalVector.z)
                            FinalVector.z = OriginalVector.z + GradientZ(OriginalVector.x, OriginalVector.y, OriginalVector.z)
                            OriginalVector.x = FinalVector.x
                            OriginalVector.y = FinalVector.y
                            OriginalVector.z = FinalVector.z
                        Loop

                    End If
                Next
            Next
        Next
        FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(Thickness, FormMain.AnatomyVolume), "Step (CSF-GM Edge) Thickness", True, True)


    End Sub
#Region "Extract Maps & Matrices"
    Public Function ReturnWMMap() As MatrixDataStructures.Matrix3DInt16
        Dim temp As New MatrixDataStructures.Matrix3DInt16(ClassifierValues.Data.x, ClassifierValues.Data.y, ClassifierValues.Data.z)
        For i As Integer = 0 To ClassifierValues.GM.x - 1
            For j As Integer = 0 To ClassifierValues.GM.y - 1
                For k As Integer = 0 To ClassifierValues.GM.z - 1
                    If ClassifierValues.GM(i, j, k) > 0 Then temp(i, j, k) = 240
                Next
            Next
        Next
        Return temp
    End Function
    Public Function ReturnGMMap() As MatrixDataStructures.Matrix3DInt16
        Dim temp As New MatrixDataStructures.Matrix3DInt16(ClassifierValues.Data.x, ClassifierValues.Data.y, ClassifierValues.Data.z)
        For i As Integer = 0 To ClassifierValues.GM.x - 1
            For j As Integer = 0 To ClassifierValues.GM.y - 1
                For k As Integer = 0 To ClassifierValues.GM.z - 1
                    If ClassifierValues.GM(i, j, k) > 0 Then temp(i, j, k) = 128
                Next
            Next
        Next
        Return temp
    End Function
    Public Function ReturnCSFMap() As MatrixDataStructures.Matrix3DInt16
        Dim temp As New MatrixDataStructures.Matrix3DInt16(ClassifierValues.Data.x, ClassifierValues.Data.y, ClassifierValues.Data.z)
        For i As Integer = 0 To ClassifierValues.GM.x - 1
            For j As Integer = 0 To ClassifierValues.GM.y - 1
                For k As Integer = 0 To ClassifierValues.GM.z - 1
                    If ClassifierValues.CSF(i, j, k) > 0 Then temp(i, j, k) = 5
                Next
            Next
        Next
        Return temp
    End Function
    Public Function ReturnGMData() As MatrixDataStructures.Matrix3DInt16
        Dim temp As New MatrixDataStructures.Matrix3DInt16(ClassifierValues.Data.x, ClassifierValues.Data.y, ClassifierValues.Data.z)
        For i As Integer = 0 To ClassifierValues.GM.x - 1
            For j As Integer = 0 To ClassifierValues.GM.y - 1
                For k As Integer = 0 To ClassifierValues.GM.z - 1
                    If ClassifierValues.GM(i, j, k) > 0 Then temp(i, j, k) = ClassifierValues.Data(i, j, k)
                Next
            Next
        Next
        Return temp
    End Function
    Public Function ReturnWMData() As MatrixDataStructures.Matrix3DInt16
        Dim temp As New MatrixDataStructures.Matrix3DInt16(ClassifierValues.Data.x, ClassifierValues.Data.y, ClassifierValues.Data.z)
        For i As Integer = 0 To ClassifierValues.GM.x - 1
            For j As Integer = 0 To ClassifierValues.GM.y - 1
                For k As Integer = 0 To ClassifierValues.GM.z - 1
                    If ClassifierValues.WM(i, j, k) > 0 Then temp(i, j, k) = ClassifierValues.Data(i, j, k)
                Next
            Next
        Next
        Return temp
    End Function
    Public Function ReturnCSFData() As MatrixDataStructures.Matrix3DInt16
        Dim temp As New MatrixDataStructures.Matrix3DInt16(ClassifierValues.Data.x, ClassifierValues.Data.y, ClassifierValues.Data.z)
        For i As Integer = 0 To ClassifierValues.GM.x - 1
            For j As Integer = 0 To ClassifierValues.GM.y - 1
                For k As Integer = 0 To ClassifierValues.GM.z - 1
                    If ClassifierValues.CSF(i, j, k) > 0 Then temp(i, j, k) = ClassifierValues.Data(i, j, k)
                Next
            Next
        Next
        Return temp
    End Function

    Public Sub PromptSaveWM()
        Dim FileSaveObject As New MatrixDataStructures.SerializeMatrixtoXML
        Dim TempVolume As New MatrixDataStructures.SingleVolume
        TempVolume.Center = FormMain.AnatomyVolume.Center
        TempVolume.FOV = FormMain.AnatomyVolume.FOV
        TempVolume.HeaderFileName = FormMain.AnatomyVolume.HeaderFileName
        If MsgBox("Save WM Probability Map", MsgBoxStyle.YesNo, "Probability Map") = MsgBoxResult.Yes Then
            TempVolume.Matrix = Me.ReturnWMMap
            FileSaveObject.SaveAs(TempVolume)
        End If
        If MsgBox("Save WM Data Map", MsgBoxStyle.YesNo, "Data Map") = MsgBoxResult.Yes Then
            TempVolume.Matrix = Me.ReturnWMData
            FileSaveObject.SaveAs(TempVolume)
        End If
    End Sub
    Public Sub PromptSaveGM()
        Dim FileSaveObject As New MatrixDataStructures.SerializeMatrixtoXML
        Dim TempVolume As New MatrixDataStructures.SingleVolume
        TempVolume.Center = FormMain.AnatomyVolume.Center
        TempVolume.FOV = FormMain.AnatomyVolume.FOV
        TempVolume.HeaderFileName = FormMain.AnatomyVolume.HeaderFileName
        If MsgBox("Save GM Probability Map", MsgBoxStyle.YesNo, "Probability Map") = MsgBoxResult.Yes Then
            TempVolume.Matrix = Me.ReturnGMMap
            FileSaveObject.SaveAs(TempVolume)
        End If
        If MsgBox("Save GM Data Map", MsgBoxStyle.YesNo, "Data Map") = MsgBoxResult.Yes Then
            TempVolume.Matrix = Me.ReturnGMData
            FileSaveObject.SaveAs(TempVolume)
        End If
    End Sub
    Public Sub PromptSaveCSF()
        Dim FileSaveObject As New MatrixDataStructures.SerializeMatrixtoXML
        Dim TempVolume As New MatrixDataStructures.SingleVolume
        TempVolume.Center = FormMain.AnatomyVolume.Center
        TempVolume.FOV = FormMain.AnatomyVolume.FOV
        TempVolume.HeaderFileName = FormMain.AnatomyVolume.HeaderFileName
        If MsgBox("Save CSF Probability Map", MsgBoxStyle.YesNo, "Probability Map") = MsgBoxResult.Yes Then
            TempVolume.Matrix = Me.ReturnCSFMap
            FileSaveObject.SaveAs(TempVolume)
        End If
        If MsgBox("Save CSF Data Map", MsgBoxStyle.YesNo, "Data Map") = MsgBoxResult.Yes Then
            TempVolume.Matrix = Me.ReturnCSFData
            FileSaveObject.SaveAs(TempVolume)
        End If
    End Sub

    Public Sub OpenGM()
        FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(Me.ReturnGMData, FormMain.AnatomyVolume), "GM Data", True, True)
        FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(Me.ReturnGMMap, FormMain.AnatomyVolume), "GM Map", True, False)
    End Sub
    Public Sub OpenWM()
        FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(Me.ReturnWMData, FormMain.AnatomyVolume), "WM Data", True, False)
        FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(Me.ReturnWMMap, FormMain.AnatomyVolume), "WM Map", True, False)
    End Sub
    Public Sub OpenCSF()
        FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(Me.ReturnCSFData, FormMain.AnatomyVolume), "CSF Data", True, False)
        FormMain.AnatomyDisplay.OpenFile(New MatrixDataStructures.SingleVolume(Me.ReturnCSFMap, FormMain.AnatomyVolume), "CSF Map", True, False)
    End Sub
#End Region

End Class

#Region " Extra Filter Classes"

Public Class DragTree
    Inherits TreeView
    Protected _leftDown As Boolean
    Protected _leftPoint As System.Drawing.Point
    Protected _leftNode As TreeNode

    Sub New()
        _leftDown = False
    End Sub

    Protected Overrides Sub OnMouseDown(ByVal e As MouseEventArgs)
        ' Only interested in the left button
        If (e.Button = Windows.Forms.MouseButtons.Left) Then
            Dim n As TreeNode = Me.GetNodeAt(New Point(e.X, e.Y))

            ' Are we selecting a valid node?
            If Not (n Is Nothing) Then
                ' Might be start of a drag, so remember details
                _leftNode = n
                _leftDown = True
                _leftPoint = New Point(e.X, e.Y)

                ' Must capture the mouse
                Me.Capture = True
                Me.Focus()
            End If
        End If

        MyBase.OnMouseDown(e)
    End Sub

    Protected Overrides Sub OnMouseMove(ByVal e As MouseEventArgs)
        ' Are we monitoring for a drag operation?
        If (_leftDown) Then
            Dim dragRect As Rectangle = New Rectangle(_leftPoint, New Size(0, 0))

            ' Create rectangle for drag start
            dragRect.Inflate(SystemInformation.DoubleClickSize)

            ' Has mouse been dragged outside of rectangle?
            If (Not dragRect.Contains(New Point(e.X, e.Y))) Then
                ' Create an object the TabbedGroups control understands
                Dim dp As TabbedGroups.DragProvider = New TabbedGroups.DragProvider()

                ' Box the node name as the parameter for passing across
                dp.Tag = CType(_leftNode.Text, Object)

                ' Must start a drag operation
                DoDragDrop(dp, DragDropEffects.Copy)

                ' Cancel any further drag events until mouse is pressed again
                _leftDown = False
                _leftNode = Nothing
            End If
        End If

        MyBase.OnMouseMove(e)
    End Sub

    Protected Overrides Sub WndProc(ByRef m As Message)
        Select Case m.Msg
            Case CType(Crownwood.Magic.Win32.Msgs.WM_LBUTTONUP, Int32)
                ' Remembering drag info?
                If (_leftDown) Then
                    ' Cancel any drag attempt
                    _leftDown = False
                    _leftNode = Nothing
                End If
        End Select

        MyBase.WndProc(m)
    End Sub
End Class

Public Class SobelEdgeDetect
    Public Sub New()
    End Sub
    Public Function Detect(ByVal Data As MatrixDataStructures.Matrix3DInt16) As MatrixDataStructures.Matrix3DInt16
        Dim temp As New MatrixDataStructures.Matrix3DInt16(Data.x, Data.y, Data.z)
        Dim a As New MorphologyFilters
        Dim Gx, Gy, G As Integer
        For k As Integer = 1 To Data.z - 2
            'work on transverse slices
            For i As Integer = 1 To Data.x - 2
                For j As Integer = 1 To Data.y - 2
                    'skip boundary edges
                    Gx = -1 * Data(i - 1, j - 1, k + 0) + 1 * Data(i + 1, j - 1, k + 0) + _
                         -2 * Data(i - 1, j + 0, k + 0) + 2 * Data(i + 1, j + 0, k + 0) + _
                         -1 * Data(i - 1, j + 1, k + 0) + 1 * Data(i + 1, j + 1, k + 0)
                    Gy = 1 * Data(i - 1, j - 1, k + 0) + 2 * Data(i, j - 1, k + 0) + 1 * Data(i + 1, j - 1, k + 0) + _
                        -1 * Data(i - 1, j + 1, k + 0) - 2 * Data(i, j + 1, k + 0) - 1 * Data(i + 1, j + 1, k + 0)
                    G = Math.Abs(Gx) + Math.Abs(Gy)
                    If G < 0 Then G = 0
                    temp(i, j, k) = CShort(G)
                Next
            Next
        Next
        Return temp
    End Function
    Public Function DetectEgge3Doptimized(ByVal Data As MatrixDataStructures.Matrix3DInt16) As MatrixDataStructures.Matrix3DInt16
        Dim temp As New MatrixDataStructures.Matrix3DInt16(Data.x, Data.y, Data.z)
        Dim a As New MorphologyFilters
        Dim Gx, Gy, G As Integer
        For k As Integer = 1 To Data.z - 2
            'work on transverse slices
            For i As Integer = 1 To Data.x - 2
                For j As Integer = 1 To Data.y - 2
                    'skip boundary edges
                    Gx = -1 * Data(i - 1, j - 1, k + 0) + 1 * Data(i + 1, j - 1, k + 0) + _
                         -2 * Data(i - 1, j + 0, k + 0) + 2 * Data(i + 1, j + 0, k + 0) + _
                         -1 * Data(i - 1, j + 1, k + 0) + 1 * Data(i + 1, j + 1, k + 0)
                    Gy = 1 * Data(i - 1, j - 1, k + 0) + 2 * Data(i, j - 1, k + 0) + 1 * Data(i + 1, j - 1, k + 0) + _
                        -1 * Data(i - 1, j + 1, k + 0) - 2 * Data(i, j + 1, k + 0) - 1 * Data(i + 1, j + 1, k + 0)
                    G = Math.Abs(Gx) + Math.Abs(Gy)
                    If G < 0 Then G = 0
                    temp(i, j, k) = CShort(G)
                Next
            Next
        Next
        Return temp
    End Function
End Class

Public Class MatrixAntiAliasing
    Public Sub AntiAlias(ByVal from_pic As Bitmap, ByRef to_pic As Bitmap)
        ' Get the source Bitmap.
        Dim from_bm As New Bitmap(from_pic)
        ' Make the destination Bitmap.
        Dim to_bm As New Bitmap(from_bm.Width, from_bm.Height)

        ' Copy the image.
        Dim gr As Graphics = Graphics.FromImage(to_bm)
        gr.SmoothingMode = Drawing2D.SmoothingMode.AntiAlias
        gr.PixelOffsetMode = Drawing2D.PixelOffsetMode.HighQuality
        gr.DrawImage(from_bm, 0, 0)
        ' Display the result.
        to_pic = to_bm
    End Sub
End Class

Public Class MorphologyFilters
    Public Sub New()
    End Sub
    Public Function Dilate(ByVal BinaryData As MatrixDataStructures.Matrix3DInt16, Optional ByVal use6side As Boolean = True) As MatrixDataStructures.Matrix3DInt16
        Dim Temp As New MatrixDataStructures.Matrix3DInt16(BinaryData.x, BinaryData.y, BinaryData.z)
        If use6side = False Then
            For k As Integer = 1 To BinaryData.z - 2
                'work on transverse slices
                For i As Integer = 1 To BinaryData.x - 2
                    For j As Integer = 1 To BinaryData.y - 2
                        'skip boundary edges
                        Temp(i, j, k) = BinaryData(i, j, k)
                        For ki As Integer = -1 To 1
                            For kj As Integer = -1 To 1
                                For kk As Integer = -1 To 1
                                    Temp(i, j, k) = Math.Max(Temp(i, j, k), BinaryData(i + ki, j + kj, k + kk))
                                Next
                            Next
                        Next
                    Next
                Next
            Next
        Else
            For k As Integer = 1 To BinaryData.z - 2
                'work on transverse slices
                For i As Integer = 1 To BinaryData.x - 2
                    For j As Integer = 1 To BinaryData.y - 2
                        'skip boundary edges
                        Temp(i, j, k) = BinaryData(i, j, k)
                        Temp(i, j, k) = Math.Max(Temp(i, j, k), BinaryData(i - 1, j, k))
                        Temp(i, j, k) = Math.Max(Temp(i, j, k), BinaryData(i + 1, j, k))
                        Temp(i, j, k) = Math.Max(Temp(i, j, k), BinaryData(i, j - 1, k))
                        Temp(i, j, k) = Math.Max(Temp(i, j, k), BinaryData(i, j + 1, k))
                        Temp(i, j, k) = Math.Max(Temp(i, j, k), BinaryData(i, j, k - 1))
                        Temp(i, j, k) = Math.Max(Temp(i, j, k), BinaryData(i, j, k + 1))
                    Next
                Next
            Next
        End If
        Return Temp
    End Function
    Public Function Erode(ByVal BinaryData As MatrixDataStructures.Matrix3DInt16, Optional ByVal use6side As Boolean = True) As MatrixDataStructures.Matrix3DInt16
        Dim Temp As New MatrixDataStructures.Matrix3DInt16(BinaryData.x, BinaryData.y, BinaryData.z)
        If use6side = False Then
            For k As Integer = 1 To BinaryData.z - 2
                'work on transverse slices
                For i As Integer = 1 To BinaryData.x - 2
                    For j As Integer = 1 To BinaryData.y - 2
                        'skip boundary edges
                        Temp(i, j, k) = BinaryData(i, j, k)
                        For ki As Integer = -1 To 1
                            For kj As Integer = -1 To 1
                                For kk As Integer = -1 To 1
                                    Temp(i, j, k) = Math.Min(Temp(i, j, k), BinaryData(i + ki, j + kj, k + kk))
                                Next
                            Next
                        Next
                    Next
                Next
            Next
        Else
            For k As Integer = 1 To BinaryData.z - 2
                'work on transverse slices
                For i As Integer = 1 To BinaryData.x - 2
                    For j As Integer = 1 To BinaryData.y - 2
                        'skip boundary edges
                        Temp(i, j, k) = BinaryData(i, j, k)
                        Temp(i, j, k) = Math.Min(Temp(i, j, k), BinaryData(i - 1, j, k))
                        Temp(i, j, k) = Math.Min(Temp(i, j, k), BinaryData(i + 1, j, k))
                        Temp(i, j, k) = Math.Min(Temp(i, j, k), BinaryData(i, j - 1, k))
                        Temp(i, j, k) = Math.Min(Temp(i, j, k), BinaryData(i, j + 1, k))
                        Temp(i, j, k) = Math.Min(Temp(i, j, k), BinaryData(i, j, k - 1))
                        Temp(i, j, k) = Math.Min(Temp(i, j, k), BinaryData(i, j, k + 1))
                    Next
                Next
            Next
        End If
        Return Temp
    End Function
    Public Function OldErode(ByVal BinaryData As MatrixDataStructures.Matrix3DInt16, Optional ByVal size As Integer = 3, Optional ByVal keepconnections As Boolean = False) As MatrixDataStructures.Matrix3DInt16
        Dim temp As New MatrixDataStructures.Matrix3DInt16(BinaryData.x, BinaryData.y, BinaryData.z)
        Dim count As Integer = 0
        Dim min As Integer = 0
        For k As Integer = 0 To BinaryData.z - size - 1
            'work on transverse slices
            For i As Integer = 0 To BinaryData.x - size - 1
                For j As Integer = 0 To BinaryData.y - size - 1
                    'skip boundary edges
                    temp(i, j, k) = BinaryData(i, j, k)
                    min = 0
                    count = 0
                    For ki As Integer = 0 To size - 1
                        For kj As Integer = 0 To size - 1
                            For kk As Integer = 0 To size - 1
                                If keepconnections = False Then
                                    temp(i, j, k) = Math.Min(temp(i, j, k), BinaryData(i + ki, j + kj, k + kk))
                                Else
                                    If BinaryData(i + ki, j + kj, k + kk) >= 1 Then count += 1
                                    min = Math.Min(min, BinaryData(i + ki, j + kj, k + kk))
                                End If
                            Next
                        Next
                    Next
                    If keepconnections = True And count > 5 Then
                        temp(i, j, k) = CShort(min)
                    End If
                Next
            Next
        Next
        Return temp
    End Function
    Public Function AreaOpen(ByVal BinaryData As MatrixDataStructures.Matrix3DInt16, ByVal count As Integer) As MatrixDataStructures.Matrix3DInt16
        Dim temp As New MatrixDataStructures.Matrix3DInt16(BinaryData.x, BinaryData.y, BinaryData.z)
        Dim counter As Integer
        For k As Integer = 1 To BinaryData.z - 2
            'work on transverse slices
            For i As Integer = 1 To BinaryData.x - 2
                For j As Integer = 1 To BinaryData.y - 2
                    'skip boundary edges
                    counter = 0
                    For ki As Integer = -1 To 1
                        For kj As Integer = -1 To 1
                            For kk As Integer = -1 To 1
                                If BinaryData(i + ki, j + kj, k + kk) > 0 Then counter += 1
                            Next
                        Next
                    Next
                    If counter >= count Then temp(i, j, k) = BinaryData(i, j, k)
                Next
            Next
        Next
        Return temp
    End Function
    Public Function Opening(ByVal BinaryData As MatrixDataStructures.Matrix3DInt16) As MatrixDataStructures.Matrix3DInt16
        Dim Temp As MatrixDataStructures.Matrix3DInt16
        Temp = Erode(BinaryData)
        Temp = Dilate(Temp)
        Return Temp
    End Function
    Public Function Closing(ByVal BinaryData As MatrixDataStructures.Matrix3DInt16) As MatrixDataStructures.Matrix3DInt16
        Dim Temp As MatrixDataStructures.Matrix3DInt16
        Temp = Dilate(BinaryData)
        Temp = Erode(Temp)
        Return Temp
    End Function
End Class

Public Class _3DFilters
    Public Sub New()
    End Sub
    Public Function GaussianSmoothen(ByVal Data As MatrixDataStructures.Matrix3DInt16) As MatrixDataStructures.Matrix3DInt16
        Dim Temp As New MatrixDataStructures.Matrix3DInt16(Data.x, Data.y, Data.z)
        Dim value As Double
        For k As Integer = 1 To Data.z - 2
            'work on transverse slices
            For i As Integer = 1 To Data.x - 2
                For j As Integer = 1 To Data.y - 2
                    'skip boundary edges
                    value = (Data(i - 1, j - 1, k - 1) * 0) + (Data(i + 0, j - 1, k - 1) * 1) + (Data(i + 1, j - 1, k - 1) * 0) + _
                            (Data(i - 1, j + 0, k - 1) * 1) + (Data(i + 0, j + 0, k - 1) * 3) + (Data(i + 1, j - 0, k - 1) * 1) + _
                            (Data(i - 1, j + 1, k - 1) * 0) + (Data(i + 0, j + 1, k - 1) * 1) + (Data(i + 1, j + 1, k - 1) * 0) + _
                            (Data(i - 1, j - 1, k - 0) * 1) + (Data(i + 0, j - 1, k - 0) * 3) + (Data(i + 1, j - 1, k - 0) * 1) + _
                            (Data(i - 1, j + 0, k - 0) * 3) + (Data(i + 0, j + 0, k - 0) * 9) + (Data(i + 1, j - 0, k - 0) * 3) + _
                            (Data(i - 1, j + 1, k - 0) * 1) + (Data(i + 0, j + 1, k - 0) * 3) + (Data(i + 1, j + 1, k - 0) * 1) + _
                            (Data(i - 1, j - 1, k + 1) * 0) + (Data(i + 0, j - 1, k + 1) * 1) + (Data(i + 1, j - 1, k + 1) * 0) + _
                            (Data(i - 1, j + 0, k + 1) * 1) + (Data(i + 0, j + 0, k + 1) * 3) + (Data(i + 1, j - 0, k + 1) * 1) + _
                            (Data(i - 1, j + 1, k + 1) * 0) + (Data(i + 0, j + 1, k + 1) * 1) + (Data(i + 1, j + 1, k + 1) * 0)
                    value /= 39
                    Temp(i, j, k) = CShort(value)
                Next
            Next
        Next
        Return Temp
    End Function
    Public Function Sobel3D(ByVal Data As MatrixDataStructures.Matrix3DInt16) As MatrixDataStructures.Matrix3DInt16
        Dim Temp As New MatrixDataStructures.Matrix3DInt16(Data.x, Data.y, Data.z)
        Dim Gx, Gy, Gz As Double
        For k As Integer = 1 To Data.z - 2
            'work on transverse slices
            For i As Integer = 1 To Data.x - 2
                For j As Integer = 1 To Data.y - 2

                    Gz = (Data(i - 1, j - 1, k - 1) * 1) + (Data(i + 0, j - 1, k - 1) * 2) + (Data(i + 1, j - 1, k - 1) * 1) + _
                         (Data(i - 1, j + 0, k - 1) * 2) + (Data(i + 0, j + 0, k - 1) * 4) + (Data(i + 1, j - 0, k - 1) * 2) + _
                         (Data(i - 1, j + 1, k - 1) * 1) + (Data(i + 0, j + 1, k - 1) * 2) + (Data(i + 1, j + 1, k - 1) * 1) + _
                         (Data(i - 1, j - 1, k - 0) * 0) + (Data(i + 0, j - 1, k - 0) * 0) + (Data(i + 1, j - 1, k - 0) * 0) + _
                         (Data(i - 1, j + 0, k - 0) * 0) + (Data(i + 0, j + 0, k - 0) * 0) + (Data(i + 1, j - 0, k - 0) * 0) + _
                         (Data(i - 1, j + 1, k - 0) * 0) + (Data(i + 0, j + 1, k - 0) * 0) + (Data(i + 1, j + 1, k - 0) * 0) + _
                         (Data(i - 1, j - 1, k + 1) * -1) + (Data(i + 0, j - 1, k + 1) * -2) + (Data(i + 1, j - 1, k + 1) * -1) + _
                         (Data(i - 1, j + 0, k + 1) * -2) + (Data(i + 0, j + 0, k + 1) * -4) + (Data(i + 1, j - 0, k + 1) * -2) + _
                         (Data(i - 1, j + 1, k + 1) * -1) + (Data(i + 0, j + 1, k + 1) * -2) + (Data(i + 1, j + 1, k + 1) * -1)

                    Gy = (Data(i - 1, j - 1, k - 1) * 1) + (Data(i + 0, j - 1, k - 1) * 0) + (Data(i + 1, j - 1, k - 1) * -1) + _
                         (Data(i - 1, j + 0, k - 1) * 2) + (Data(i + 0, j + 0, k - 1) * 0) + (Data(i + 1, j - 0, k - 1) * -2) + _
                         (Data(i - 1, j + 1, k - 1) * 1) + (Data(i + 0, j + 1, k - 1) * 0) + (Data(i + 1, j + 1, k - 1) * -1) + _
                         (Data(i - 1, j - 1, k - 0) * 2) + (Data(i + 0, j - 1, k - 0) * 0) + (Data(i + 1, j - 1, k - 0) * -2) + _
                         (Data(i - 1, j + 0, k - 0) * 4) + (Data(i + 0, j + 0, k - 0) * 0) + (Data(i + 1, j - 0, k - 0) * -4) + _
                         (Data(i - 1, j + 1, k - 0) * 2) + (Data(i + 0, j + 1, k - 0) * 0) + (Data(i + 1, j + 1, k - 0) * -4) + _
                         (Data(i - 1, j - 1, k + 1) * 1) + (Data(i + 0, j - 1, k + 1) * 0) + (Data(i + 1, j - 1, k + 1) * -1) + _
                         (Data(i - 1, j + 0, k + 1) * 2) + (Data(i + 0, j + 0, k + 1) * 0) + (Data(i + 1, j - 0, k + 1) * -2) + _
                         (Data(i - 1, j + 1, k + 1) * 1) + (Data(i + 0, j + 1, k + 1) * 0) + (Data(i + 1, j + 1, k + 1) * -1)

                    Gz = (Data(i - 1, j - 1, k - 1) * 1) + (Data(i + 0, j - 1, k - 1) * 2) + (Data(i + 1, j - 1, k - 1) * 1) + _
                         (Data(i - 1, j + 0, k - 1) * 0) + (Data(i + 0, j + 0, k - 1) * 0) + (Data(i + 1, j - 0, k - 1) * 0) + _
                         (Data(i - 1, j + 1, k - 1) * -1) + (Data(i + 0, j + 1, k - 1) * -2) + (Data(i + 1, j + 1, k - 1) * -1) + _
                         (Data(i - 1, j - 1, k - 0) * 2) + (Data(i + 0, j - 1, k - 0) * 4) + (Data(i + 1, j - 1, k - 0) * 2) + _
                         (Data(i - 1, j + 0, k - 0) * 0) + (Data(i + 0, j + 0, k - 0) * 0) + (Data(i + 1, j - 0, k - 0) * 0) + _
                         (Data(i - 1, j + 1, k - 0) * -2) + (Data(i + 0, j + 1, k - 0) * -4) + (Data(i + 1, j + 1, k - 0) * -2) + _
                         (Data(i - 1, j - 1, k + 1) * 1) + (Data(i + 0, j - 1, k + 1) * 2) + (Data(i + 1, j - 1, k + 1) * -1) + _
                         (Data(i - 1, j + 0, k + 1) * 0) + (Data(i + 0, j + 0, k + 1) * 0) + (Data(i + 1, j - 0, k + 1) * 0) + _
                         (Data(i - 1, j + 1, k + 1) * -1) + (Data(i + 0, j + 1, k + 1) * -2) + (Data(i + 1, j + 1, k + 1) * -1)
                    Temp(i, j, k) = CShort(Math.Abs(Gx) + Math.Abs(Gy) + Math.Abs(Gz))
                Next
            Next
        Next
        Return Temp
    End Function
    Public Function Sobel3DNonMaximumSupression(ByVal Data As MatrixDataStructures.Matrix3DInt16) As MatrixDataStructures.Matrix3DInt16
        Dim Temp As New MatrixDataStructures.Matrix3DInt16(Data.x, Data.y, Data.z)
        Dim Gx, Gy, Gz As New MatrixDataStructures.Matrix3DInt16(Data.x, Data.y, Data.z)
        For k As Integer = 1 To Data.z - 2
            For i As Integer = 1 To Data.x - 2
                For j As Integer = 1 To Data.y - 2
                    Gz(i, j, k) = CShort((Data(i - 1, j - 1, k - 1) * 1) + (Data(i + 0, j - 1, k - 1) * 2) + (Data(i + 1, j - 1, k - 1) * 1) + _
                         (Data(i - 1, j + 0, k - 1) * 2) + (Data(i + 0, j + 0, k - 1) * 4) + (Data(i + 1, j - 0, k - 1) * 2) + _
                         (Data(i - 1, j + 1, k - 1) * 1) + (Data(i + 0, j + 1, k - 1) * 2) + (Data(i + 1, j + 1, k - 1) * 1) + _
                         (Data(i - 1, j - 1, k - 0) * 0) + (Data(i + 0, j - 1, k - 0) * 0) + (Data(i + 1, j - 1, k - 0) * 0) + _
                         (Data(i - 1, j + 0, k - 0) * 0) + (Data(i + 0, j + 0, k - 0) * 0) + (Data(i + 1, j - 0, k - 0) * 0) + _
                         (Data(i - 1, j + 1, k - 0) * 0) + (Data(i + 0, j + 1, k - 0) * 0) + (Data(i + 1, j + 1, k - 0) * 0) + _
                         (Data(i - 1, j - 1, k + 1) * -1) + (Data(i + 0, j - 1, k + 1) * -2) + (Data(i + 1, j - 1, k + 1) * -1) + _
                         (Data(i - 1, j + 0, k + 1) * -2) + (Data(i + 0, j + 0, k + 1) * -4) + (Data(i + 1, j - 0, k + 1) * -2) + _
                         (Data(i - 1, j + 1, k + 1) * -1) + (Data(i + 0, j + 1, k + 1) * -2) + (Data(i + 1, j + 1, k + 1) * -1))

                    Gy(i, j, k) = CShort((Data(i - 1, j - 1, k - 1) * 1) + (Data(i + 0, j - 1, k - 1) * 0) + (Data(i + 1, j - 1, k - 1) * -1) + _
                         (Data(i - 1, j + 0, k - 1) * 2) + (Data(i + 0, j + 0, k - 1) * 0) + (Data(i + 1, j - 0, k - 1) * -2) + _
                         (Data(i - 1, j + 1, k - 1) * 1) + (Data(i + 0, j + 1, k - 1) * 0) + (Data(i + 1, j + 1, k - 1) * -1) + _
                         (Data(i - 1, j - 1, k - 0) * 2) + (Data(i + 0, j - 1, k - 0) * 0) + (Data(i + 1, j - 1, k - 0) * -2) + _
                         (Data(i - 1, j + 0, k - 0) * 4) + (Data(i + 0, j + 0, k - 0) * 0) + (Data(i + 1, j - 0, k - 0) * -4) + _
                         (Data(i - 1, j + 1, k - 0) * 2) + (Data(i + 0, j + 1, k - 0) * 0) + (Data(i + 1, j + 1, k - 0) * -4) + _
                         (Data(i - 1, j - 1, k + 1) * 1) + (Data(i + 0, j - 1, k + 1) * 0) + (Data(i + 1, j - 1, k + 1) * -1) + _
                         (Data(i - 1, j + 0, k + 1) * 2) + (Data(i + 0, j + 0, k + 1) * 0) + (Data(i + 1, j - 0, k + 1) * -2) + _
                         (Data(i - 1, j + 1, k + 1) * 1) + (Data(i + 0, j + 1, k + 1) * 0) + (Data(i + 1, j + 1, k + 1) * -1))

                    Gx(i, j, k) = CShort((Data(i - 1, j - 1, k - 1) * 1) + (Data(i + 0, j - 1, k - 1) * 2) + (Data(i + 1, j - 1, k - 1) * 1) + _
                         (Data(i - 1, j + 0, k - 1) * 0) + (Data(i + 0, j + 0, k - 1) * 0) + (Data(i + 1, j - 0, k - 1) * 0) + _
                         (Data(i - 1, j + 1, k - 1) * -1) + (Data(i + 0, j + 1, k - 1) * -2) + (Data(i + 1, j + 1, k - 1) * -1) + _
                         (Data(i - 1, j - 1, k - 0) * 2) + (Data(i + 0, j - 1, k - 0) * 4) + (Data(i + 1, j - 1, k - 0) * 2) + _
                         (Data(i - 1, j + 0, k - 0) * 0) + (Data(i + 0, j + 0, k - 0) * 0) + (Data(i + 1, j - 0, k - 0) * 0) + _
                         (Data(i - 1, j + 1, k - 0) * -2) + (Data(i + 0, j + 1, k - 0) * -4) + (Data(i + 1, j + 1, k - 0) * -2) + _
                         (Data(i - 1, j - 1, k + 1) * 1) + (Data(i + 0, j - 1, k + 1) * 2) + (Data(i + 1, j - 1, k + 1) * -1) + _
                         (Data(i - 1, j + 0, k + 1) * 0) + (Data(i + 0, j + 0, k + 1) * 0) + (Data(i + 1, j - 0, k + 1) * 0) + _
                         (Data(i - 1, j + 1, k + 1) * -1) + (Data(i + 0, j + 1, k + 1) * -2) + (Data(i + 1, j + 1, k + 1) * -1))
                    Temp(i, j, k) = CShort(Math.Sqrt((Gx(i, j, k) ^ 2) + (Gy(i, j, k) ^ 2) + (Gz(i, j, k) ^ 2)))
                Next
            Next
        Next

        For k As Integer = 1 To Data.z - 2
            For i As Integer = 1 To Data.x - 2
                For j As Integer = 1 To Data.y - 2

                    If Gx(i, j, k) = Math.Max(Gx(i, j, k), Math.Max(Gy(i, j, k), Gz(i, j, k))) Then
                        If Gx(i, j, k) < Gx(i + 1, j, k) Or Gx(i, j, k) < Gx(i - 1, j, k) Then Temp(i, j, k) = 0
                    ElseIf Gy(i, j, k) = Math.Max(Gx(i, j, k), Math.Max(Gy(i, j, k), Gz(i, j, k))) Then
                        If Gx(i, j, k) < Gx(i, j + 1, k) Or Gx(i, j, k) < Gx(i, j - 1, k) Then Temp(i, j, k) = 0
                    Else
                        Gz(i, j, k) = Math.Max(Gx(i, j, k), Math.Max(Gy(i, j, k), Gz(i, j, k)))
                        If Gx(i, j, k) < Gx(i, j, k + 1) Or Gx(i, j, k) < Gx(i, j, k - 1) Then Temp(i, j, k) = 0
                    End If
                Next
            Next
        Next
        Return Temp
    End Function
End Class

Public Class Skeletonization
    Public Function Skeletonize(ByVal data As MatrixDataStructures.Matrix3DInt16, Optional ByVal threshold As Double = 0.2) As MatrixDataStructures.Matrix3DInt16
        Dim temp As New MatrixDataStructures.Matrix3DInt16(data.x, data.y, data.z)
        For i As Integer = 0 To data.x - 1
            For j As Integer = 0 To data.y - 1
                For k As Integer = 0 To data.z - 1
                    temp(i, j, k) = data(i, j, k)
                Next
            Next
        Next
        Dim tempm As twoDMatrixLibrary.Matrix
        For i As Integer = 0 To data.x - 1
            tempm = New twoDMatrixLibrary.Matrix(data.y, data.z)
            For j As Integer = 0 To data.y - 1
                For k As Integer = 0 To data.z - 1
                    tempm(j, k) = temp(i, j, k)
                Next
            Next
            Me.ComputeDistanceTransform2D(tempm)
            For j As Integer = 0 To data.y - 1
                For k As Integer = 0 To data.z - 1
                    temp(i, j, k) = CShort(tempm(j, k))
                Next
            Next
        Next
        Dim max As Double = temp.Maximum * threshold
        For i As Integer = 0 To data.x - 1
            For j As Integer = 0 To data.y - 1
                For k As Integer = 0 To data.z - 1
                    If temp(i, j, k) < max Then temp(i, j, k) = 0 Else temp(i, j, k) = 1
                Next
            Next
        Next
        Return temp
    End Function

    Private Sub ComputeDistanceTransform2D(ByRef data As twoDMatrixLibrary.Matrix)
        '  dt of 2d function Imports squared distance */
        Dim width As Integer = data.Rows
        Dim height As Integer = data.Columns
        Dim f(Math.Max(width, height)) As Double
        ' transform along columns
        For x As Integer = 0 To width - 1
            For y As Integer = 0 To height - 1
                f(y) = data(x, y)
            Next
            Dim d() As Double = Me.DistanceTransform1D(f, height)
            For y As Integer = 0 To height - 1
                data(x, y) = d(y)
            Next
        Next
        ' transform along rows
        For y As Integer = 0 To height - 1
            For x As Integer = 0 To width - 1 Step x + 1
                f(x) = data(x, y)
            Next
            Dim d() As Double = Me.DistanceTransform1D(f, width)
            For x As Integer = 0 To width - 1 Step x + 1
                data(x, y) = d(x)
            Next
        Next

    End Sub
    Private Function DistanceTransform1D(ByVal f As Double(), ByVal n As Integer) As Double()
        Dim d(n) As Double
        Dim v(n) As Integer
        Dim z(n + 1) As Double
        Dim k As Integer = 0

        v(0) = 0
        z(0) = -1 * 1.0E+20
        z(1) = +1 * 1.0E+20

        For q As Integer = 1 To n - 1
            Dim s As Double = ((f(q) + (q ^ 2)) - (f(v(k)) + (v(k) ^ 2))) / (2 * q - 2 * v(k))
            While (s <= z(k))
                k -= 1
                s = ((f(q) + (q ^ 2)) - (f(v(k)) + (v(k) ^ 2))) / (2 * q - 2 * v(k))
            End While
            k += 1
            v(k) = q
            z(k) = s
            z(k + 1) = +1 * 1.0E+20
        Next
        k = 0
        For q As Integer = 0 To n - 1
            While (z(k + 1) < q)
                k += 1
                d(q) = ((q - v(k)) ^ 2) + f(v(k))
            End While
        Next
        Return d

    End Function



End Class

#End Region